NC13955 QTL Analysis Script
Zach Winn
2023-05-07
Introduction
This is a reanalysis of the NC13955 population QTL analysis. This is done to produce a single, specific QTL analysis for publication.
#set working directory
setwd("C:/Users/zwinn/OneDrive/Publication/Manuscripts In Progress/FHB QTL/Final Formatted Data and Analysis")
#read in data
pheno<-read.csv("NC13955 Combined Data Final.csv")
#make sure everything is the right value
pheno[,1:7]<-lapply(pheno[,1:7], as.factor)
pheno[,8:ncol(pheno)]<-lapply(pheno[,8:ncol(pheno)], as.numeric)
#read in genetic information
geno<-read.vcf("NC13-20076xGA06493-13LE6_filt.vcf.gz", convert.chr = F)## ped stats and snps stats have been set.
## 'p' has been set.
## 'mu' and 'sigma' have been set.#convert genotype matrix to format for rqtl
geno<-format_qtlmap_geno(geno,
par_a="13955-GA06493-13LE6",
par_b="13955-NC13-20076",
rm_het=T,
rm_miss=T,
include_pars=T,
out_fmt="rqtl")$abh
#replace heterozygous calls with missing data
geno<-replace(geno, geno=="H", "-")Calculate BLUEs within and across environments
#set up lists for loop
env<-levels(pheno$ENV)
traits<-colnames(pheno)[8:12]
#create dataframe for means out
BLUEs_loc<-data.frame(GENOTYPE=unique(pheno$GENOTYPE))
#check normality
for (i in traits){
hist(pheno[[i]], main=i)
}
#run models for each environment
for (j in env){
for (i in traits){
location<-pheno[pheno$ENV==j,] #filter out each environment
print(paste("------- Analyzing trait", i, "in", j,"-------"))
if((sum(is.na(location[[i]]==T)))>200){ #check if there is data
print("Data not found, moving to the next trait")
}else
if (i=="FDK"){
print(paste("Analyzing with GMLM with Inverse Gaussian distribution"))
#Run Model
mlm<-asreml(fixed = location[,i]~1+GENOTYPE,
random = ~ (REP),
residual = ~ idv(units),
data = location,
family= asr_inverse.gaussian())
#Check significance
print(wald(mlm))
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_", j, sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}else
if(i=="DON"){
print(paste("Analyzing with GMLM with Inverse Gaussian distribution"))
#Run Model
mlm<-asreml(fixed = location[,i]~1+GENOTYPE,
random = ~ (REP),
residual = ~ idv(units),
data = location,
family= asr_inverse.gaussian())
#Check significance
print(wald(mlm))
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_", j, sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}else{
print(paste("Analyzing with MLM with Gaussian distribution"))
#Run Model
mlm<-asreml(fixed = location[,i]~1+GENOTYPE,
random = ~ (REP),
residual = ~ idv(units),
data = location)
#Check significance
print(wald(mlm))
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_", j, sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_loc<-left_join(BLUEs_loc, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}
}
}## [1] "------- Analyzing trait HD in KIN2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FHB in KIN2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:08 2023
## LogLik Sigma2 DF wall cpu
## 1 -228.014 1.0 149 13:35:08 0.0 (1 restrained)
## 2 -181.290 1.0 149 13:35:08 0.0 (1 restrained)
## 3 -133.124 1.0 149 13:35:09 0.0 (1 restrained)
## 4 -104.473 1.0 149 13:35:09 0.0 (1 restrained)
## 5 -93.401 1.0 149 13:35:09 0.0 (1 restrained)
## 6 -92.100 1.0 149 13:35:09 0.0
## 7 -92.070 1.0 149 13:35:09 0.0
## 8 -92.070 1.0 149 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 7676.9 7676.9 < 2.2e-16 ***
## GENOTYPE 184 1420.4 1420.4 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FDK in KIN2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -850.740 1.0 158 13:35:09 0.0
## 2 -757.507 1.0 158 13:35:09 0.0
## 3 -656.700 1.0 158 13:35:09 0.0
## 4 -590.195 1.0 158 13:35:09 0.0
## 5 -558.335 1.0 158 13:35:09 0.0
## 6 -551.609 1.0 158 13:35:09 0.0
## 7 -551.069 1.0 158 13:35:09 0.0
## 8 -551.063 1.0 158 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 2288.37 2288.37 < 2.2e-16 ***
## GENOTYPE 184 862.08 862.08 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait DON in KIN2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -819.956 1.0 158 13:35:09 0.0
## 2 -737.235 1.0 158 13:35:09 0.0
## 3 -648.235 1.0 158 13:35:09 0.0
## 4 -590.321 1.0 158 13:35:09 0.0
## 5 -563.424 1.0 158 13:35:09 0.0
## 6 -558.168 1.0 158 13:35:09 0.0
## 7 -557.780 1.0 158 13:35:09 0.0
## 8 -557.771 1.0 158 13:35:09 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 233.26 233.26 < 2.2e-16 ***
## GENOTYPE 184 967.91 967.91 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in KIN2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -641.368 1.0 232 13:35:09 0.0 (1 restrained)
## 2 -537.434 1.0 232 13:35:09 0.0 (1 restrained)
## 3 -427.349 1.0 232 13:35:09 0.0 (1 restrained)
## 4 -357.531 1.0 232 13:35:09 0.0 (1 restrained)
## 5 -326.764 1.0 232 13:35:09 0.0 (1 restrained)
## 6 -321.594 1.0 232 13:35:09 0.0
## 7 -321.337 1.0 232 13:35:09 0.0
## 8 -321.336 1.0 232 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1021841 1021841 < 2.2e-16 ***
## GENOTYPE 192 1327 1327 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FHB in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -882.862 1.0 231 13:35:09 0.0
## 2 -693.549 1.0 231 13:35:09 0.0
## 3 -486.358 1.0 231 13:35:09 0.0
## 4 -344.914 1.0 231 13:35:09 0.0
## 5 -271.571 1.0 231 13:35:09 0.0
## 6 -252.398 1.0 231 13:35:09 0.0
## 7 -249.968 1.0 231 13:35:09 0.0
## 8 -249.902 1.0 231 13:35:09 0.0
## 9 -249.902 1.0 231 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 352.28 352.28 < 2.2e-16 ***
## GENOTYPE 193 732.77 732.77 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait FDK in KIN2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -1169.203 1.0 233 13:35:09 0.0
## 2 -1043.001 1.0 233 13:35:09 0.0
## 3 -906.989 1.0 233 13:35:09 0.0
## 4 -818.037 1.0 233 13:35:09 0.0
## 5 -776.170 1.0 233 13:35:09 0.0
## 6 -767.534 1.0 233 13:35:09 0.0
## 7 -766.662 1.0 233 13:35:09 0.0
## 8 -766.560 1.0 233 13:35:09 0.0
## 9 -766.546 1.0 233 13:35:09 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 541.64 541.64 < 2.2e-16 ***
## GENOTYPE 193 1141.84 1141.84 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait DON in KIN2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -671.414 1.0 182 13:35:09 0.0
## 2 -600.399 1.0 182 13:35:09 0.0
## 3 -525.043 1.0 182 13:35:09 0.0
## 4 -477.866 1.0 182 13:35:09 0.0
## 5 -457.810 1.0 182 13:35:09 0.0
## 6 -454.704 1.0 182 13:35:09 0.0
## 7 -454.540 1.0 182 13:35:09 0.0
## 8 -454.534 1.0 182 13:35:09 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 190.79 190.79 < 2.2e-16 ***
## GENOTYPE 193 1325.51 1325.51 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in KIN2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -860.120 1.0 234 13:35:09 0.0 (1 restrained)
## 2 -743.455 1.0 234 13:35:09 0.0
## 3 -619.080 1.0 234 13:35:09 0.0
## 4 -538.702 1.0 234 13:35:09 0.0
## 5 -501.992 1.0 234 13:35:09 0.0
## 6 -495.193 1.0 234 13:35:09 0.0
## 7 -494.779 1.0 234 13:35:09 0.0
## 8 -494.777 1.0 234 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 86974 86974 < 2.2e-16 ***
## GENOTYPE 193 1211 1211 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait HD in LW2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -691.787 1.0 192 13:35:09 0.0 (1 restrained)
## 2 -559.886 1.0 192 13:35:09 0.0
## 3 -416.268 1.0 192 13:35:09 0.0
## 4 -319.913 1.0 192 13:35:09 0.0
## 5 -271.918 1.0 192 13:35:09 0.0
## 6 -260.652 1.0 192 13:35:09 0.0
## 7 -259.517 1.0 192 13:35:09 0.0
## 8 -259.498 1.0 192 13:35:09 0.0
## 9 -259.498 1.0 192 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 300979 300979 < 2.2e-16 ***
## GENOTYPE 186 804 804 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FHB in LW2019 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -317.568 1.0 186 13:35:09 0.0 (1 restrained)
## 2 -247.435 1.0 186 13:35:09 0.0 (1 restrained)
## 3 -173.866 1.0 186 13:35:09 0.0 (1 restrained)
## 4 -128.702 1.0 186 13:35:09 0.0 (1 restrained)
## 5 -110.014 1.0 186 13:35:09 0.0 (1 restrained)
## 6 -107.354 1.0 186 13:35:09 0.0
## 7 -107.261 1.0 186 13:35:09 0.0
## 8 -107.261 1.0 186 13:35:09 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 9558.3 9558.3 < 2.2e-16 ***
## GENOTYPE 185 1358.2 1358.2 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FDK in LW2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:09 2023
## LogLik Sigma2 DF wall cpu
## 1 -1337.441 1.0 197 13:35:10 0.0
## 2 -1153.277 1.0 197 13:35:10 0.0
## 3 -950.889 1.0 197 13:35:10 0.0
## 4 -811.154 1.0 197 13:35:10 0.0
## 5 -736.793 1.0 197 13:35:10 0.0
## 6 -715.957 1.0 197 13:35:10 0.0
## 7 -712.894 1.0 197 13:35:10 0.0
## 8 -712.779 1.0 197 13:35:10 0.0
## 9 -712.779 1.0 197 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1430.25 1430.25 < 2.2e-16 ***
## GENOTYPE 186 600.11 600.11 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait DON in LW2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -1414.045 1.0 197 13:35:10 0.0
## 2 -1216.349 1.0 197 13:35:10 0.0
## 3 -998.602 1.0 197 13:35:10 0.0
## 4 -847.293 1.0 197 13:35:10 0.0
## 5 -765.488 1.0 197 13:35:10 0.0
## 6 -741.388 1.0 197 13:35:10 0.0
## 7 -737.166 1.0 197 13:35:10 0.0
## 8 -736.727 1.0 197 13:35:10 0.0
## 9 -736.649 1.0 197 13:35:10 0.0
## 10 -736.640 1.0 197 13:35:10 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 61.23 61.23 5.107e-15 ***
## GENOTYPE 186 569.20 569.20 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in LW2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -972.328 1.0 222 13:35:10 0.0
## 2 -772.024 1.0 222 13:35:10 0.0
## 3 -552.110 1.0 222 13:35:10 0.0
## 4 -400.678 1.0 222 13:35:10 0.0
## 5 -320.581 1.0 222 13:35:10 0.0
## 6 -298.500 1.0 222 13:35:10 0.0
## 7 -295.365 1.0 222 13:35:10 0.0
## 8 -295.256 1.0 222 13:35:10 0.0
## 9 -295.256 1.0 222 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 40395 40395 < 2.2e-16 ***
## GENOTYPE 192 635 635 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FHB in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -482.340 1.0 223 13:35:10 0.0 (1 restrained)
## 2 -381.960 1.0 223 13:35:10 0.0 (1 restrained)
## 3 -275.882 1.0 223 13:35:10 0.0
## 4 -208.294 1.0 223 13:35:10 0.0
## 5 -178.416 1.0 223 13:35:10 0.0 (1 restrained)
## 6 -173.362 1.0 223 13:35:10 0.0 (1 restrained)
## 7 -173.108 1.0 223 13:35:10 0.0
## 8 -173.107 1.0 223 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 6294.6 6294.6 < 2.2e-16 ***
## GENOTYPE 192 1294.8 1294.8 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait FDK in LW2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -1079.621 1.0 218 13:35:10 0.0
## 2 -958.974 1.0 218 13:35:10 0.0
## 3 -828.869 1.0 218 13:35:10 0.0
## 4 -743.660 1.0 218 13:35:10 0.0
## 5 -703.515 1.0 218 13:35:10 0.0
## 6 -695.396 1.0 218 13:35:10 0.0
## 7 -694.784 1.0 218 13:35:10 0.0
## 8 -694.776 1.0 218 13:35:10 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 444.85 444.85 < 2.2e-16 ***
## GENOTYPE 192 1070.34 1070.34 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in LW2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -356.315 1.0 106 13:35:10 0.0
## 2 -302.485 1.0 106 13:35:10 0.0
## 3 -244.630 1.0 106 13:35:10 0.0
## 4 -207.098 1.0 106 13:35:10 0.0
## 5 -189.805 1.0 106 13:35:10 0.0
## 6 -186.524 1.0 106 13:35:10 0.0
## 7 -186.313 1.0 106 13:35:10 0.0
## 8 -186.312 1.0 106 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 438.52 438.52 < 2.2e-16 ***
## GENOTYPE 184 851.54 851.54 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in LW2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -902.787 1.0 223 13:35:10 0.0 (1 restrained)
## 2 -765.630 1.0 223 13:35:10 0.0
## 3 -618.001 1.0 223 13:35:10 0.0
## 4 -520.171 1.0 223 13:35:10 0.0
## 5 -472.821 1.0 223 13:35:10 0.0
## 6 -462.545 1.0 223 13:35:10 0.0
## 7 -461.669 1.0 223 13:35:10 0.0
## 8 -461.658 1.0 223 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 130425 130425 < 2.2e-16 ***
## GENOTYPE 191 976 976 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait HD in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FHB in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait FDK in VIR2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -1020.821 1.0 195 13:35:10 0.0
## 2 -922.179 1.0 195 13:35:10 0.0
## 3 -816.221 1.0 195 13:35:10 0.0
## 4 -747.594 1.0 195 13:35:10 0.0
## 5 -716.091 1.0 195 13:35:10 0.0
## 6 -710.172 1.0 195 13:35:10 0.0
## 7 -709.799 1.0 195 13:35:10 0.0
## 8 -709.796 1.0 195 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1028.7 1028.7 < 2.2e-16 ***
## GENOTYPE 186 1080.5 1080.5 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "------- Analyzing trait DON in VIR2019 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -576.052 1.0 196 13:35:10 0.0
## 2 -535.425 1.0 196 13:35:10 0.0
## 3 -493.851 1.0 196 13:35:10 0.0
## 4 -470.318 1.0 196 13:35:10 0.0
## 5 -462.305 1.0 196 13:35:10 0.0
## 6 -461.553 1.0 196 13:35:10 0.0
## 7 -461.504 1.0 196 13:35:10 0.0
## 8 -461.500 1.0 196 13:35:10 0.0## Warning in asreml(fixed = location[, i] ~ 1 + GENOTYPE, random = ~(REP), : Some
## components changed by more than 1% on the last iteration.
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 148.11 148.11 < 2.2e-16 ***
## GENOTYPE 186 2082.22 2082.22 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in VIR2019 -------"
## [1] "Data not found, moving to the next trait"
## [1] "------- Analyzing trait HD in VIR2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:10 2023
## LogLik Sigma2 DF wall cpu
## 1 -611.866 1.0 234 13:35:10 0.0 (1 restrained)
## 2 -495.363 1.0 234 13:35:10 0.0
## 3 -371.420 1.0 234 13:35:10 0.0
## 4 -291.397 1.0 234 13:35:10 0.0
## 5 -254.889 1.0 234 13:35:10 0.0
## 6 -248.148 1.0 234 13:35:10 0.0
## 7 -247.741 1.0 234 13:35:10 0.0
## 8 -247.738 1.0 234 13:35:10 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 2559200 2559200 < 2.2e-16 ***
## GENOTYPE 194 1219 1219 < 2.2e-16 ***
## residual (MS) 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FHB in VIR2020 -------"
## [1] "Analyzing with MLM with Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:11 2023
## LogLik Sigma2 DF wall cpu
## 1 -614.597 1.0 233 13:35:11 0.0 (1 restrained)
## 2 -492.704 1.0 233 13:35:11 0.0
## 3 -362.606 1.0 233 13:35:11 0.0
## 4 -278.034 1.0 233 13:35:11 0.0
## 5 -238.879 1.0 233 13:35:11 0.0
## 6 -231.352 1.0 233 13:35:11 0.0
## 7 -230.856 1.0 233 13:35:11 0.0
## 8 -230.853 1.0 233 13:35:11 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1047.2 1047.2 < 2.2e-16 ***
## GENOTYPE 194 1163.7 1163.7 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait FDK in VIR2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:11 2023
## LogLik Sigma2 DF wall cpu
## 1 -1116.367 1.0 228 13:35:11 0.0
## 2 -989.165 1.0 228 13:35:11 0.0
## 3 -851.930 1.0 228 13:35:11 0.0
## 4 -761.943 1.0 228 13:35:11 0.0 (1 restrained)
## 5 -719.476 1.0 228 13:35:11 0.0 (1 restrained)
## 6 -710.884 1.0 228 13:35:11 0.0 (1 restrained)
## 7 -710.252 1.0 228 13:35:11 0.0 (1 restrained)
## 8 -710.246 1.0 228 13:35:11 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 2234.4 2234.4 < 2.2e-16 ***
## GENOTYPE 194 1085.2 1085.2 < 2.2e-16 ***
## residual (MS) 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait DON in VIR2020 -------"
## [1] "Analyzing with GMLM with Inverse Gaussian distribution"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:11 2023
## LogLik Sigma2 DF wall cpu
## 1 -327.582 1.0 225 13:35:11 0.0 (1 restrained)
## 2 -272.817 1.0 225 13:35:11 0.0
## 3 -215.680 1.0 225 13:35:11 0.0
## 4 -182.480 1.0 225 13:35:11 0.0
## 5 -170.610 1.0 225 13:35:11 0.0
## 6 -169.501 1.0 225 13:35:11 0.0
## 7 -169.485 1.0 225 13:35:11 0.0
## 8 -169.485 1.0 225 13:35:11 0.0
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: location[, i][0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 552.94 552.94 < 2.2e-16 ***
## GENOTYPE 194 2020.25 2020.25 < 2.2e-16 ***
## residual (MS) 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "------- Analyzing trait PH in VIR2020 -------"
## [1] "Data not found, moving to the next trait"Multi-Environmental Models
#define traits
traits<-colnames(pheno)[8:12]
#make dataframe for predictions
BLUEs_mem<-data.frame(GENOTYPE=unique(pheno$GENOTYPE))
#run MEMLM loop
for (i in traits){
print(paste("------- Analyzing trait", i, "-------"))
siteyears<-pheno %>%
select(ENV, YEAR, LOC, all_of(i)) %>%
drop_na() %>%
distinct(YEAR)
siteyears<-as.numeric(count(siteyears))
if(((siteyears)>=2)){
print(paste(i, "has multiple years of data"))
if(i=="FDK"){
#run models
mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
random = ~ ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75,
family= asr_inverse.gaussian())
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_ME", sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}else
if(i=="DON"){
#run models
mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
random = ~ ENV+
GENOTYPE:ENV,
data = pheno,
maxit=75,
family= asr_inverse.gaussian())
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_ME", sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}else{
#run models
mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
random = ~ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75)
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_ME", sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}
}else{
print(paste(i, "has a single year of data"))
#run model
mlm<-asreml(fixed = pheno[,i]~1+GENOTYPE,
random = ~ENV+
GENOTYPE:ENV,
residual = ~units,
data = pheno,
maxit=75)
#print resid
qqnorm(mlm$residuals)
#Pull BLUEs
beta_hat<-mlm$coefficients$fixed
beta_0<-beta_hat[grep("(Intercept)", rownames(beta_hat))]
beta_1<-beta_hat[grep("GENOTYPE", rownames(beta_hat))]
BLUEs<-as.data.frame(beta_0+beta_1)
rownames(BLUEs)<-gsub("GENOTYPE_", "", rownames(beta_hat)[1:nrow(beta_hat)-1])
colnames(BLUEs)<-paste(i, "_ME", sep = "")
BLUEs<-rownames_to_column(BLUEs, var="GENOTYPE")
BLUEs_mem<-left_join(BLUEs_mem, BLUEs, by="GENOTYPE")
#Remove Junk
remove(beta_hat, beta_0, beta_1, BLUEs, mlm)
}
}## [1] "------- Analyzing trait HD -------"
## [1] "HD has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:12 2023
## LogLik Sigma2 DF wall cpu
## 1 -2381.807 1.0 1451 13:35:12 0.0 (1 restrained)
## Log-likelihood decreased to -2678.06 ; trying with reduced updates (0.81626)
## 2 -2325.628 1.0 1451 13:35:12 0.0
## 3 -2030.370 1.0 1451 13:35:12 0.0
## 4 -1798.807 1.0 1451 13:35:12 0.0
## 5 -1771.640 1.0 1451 13:35:12 0.0
## 6 -1763.747 1.0 1451 13:35:12 0.0
## 7 -1760.076 1.0 1451 13:35:12 0.0
## 8 -1758.663 1.0 1451 13:35:12 0.0
## 9 -1758.233 1.0 1451 13:35:12 0.0
## 10 -1758.157 1.0 1451 13:35:12 0.0
## 11 -1758.153 1.0 1451 13:35:12 0.0
## [1] "------- Analyzing trait FHB -------"
## [1] "FHB has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:12 2023
## LogLik Sigma2 DF wall cpu
## 1 -3215.668 1.0 1779 13:35:12 0.0
## 2 -2618.147 1.0 1779 13:35:12 0.0
## 3 -1986.339 1.0 1779 13:35:12 0.0
## 4 -1593.414 1.0 1779 13:35:12 0.0
## 5 -1426.382 1.0 1779 13:35:12 0.0
## 6 -1399.236 1.0 1779 13:35:12 0.0
## 7 -1397.875 1.0 1779 13:35:12 0.0
## 8 -1397.871 1.0 1779 13:35:12 0.0
## [1] "------- Analyzing trait FDK -------"
## [1] "FDK has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:12 2023
## LogLik Sigma2 DF wall cpu
## 1 -12793.71 1.0 2174 13:35:12 0.0
## 2 -11133.10 1.0 2174 13:35:12 0.0
## 3 -9316.39 1.0 2174 13:35:12 0.0
## 4 -8075.14 1.0 2174 13:35:12 0.0
## 5 -7421.83 1.0 2174 13:35:12 0.0
## 6 -7226.13 1.0 2174 13:35:12 0.0
## 7 -7172.62 1.0 2174 13:35:12 0.0
## 8 -7152.58 1.0 2174 13:35:12 0.0
## 9 -7143.52 1.0 2174 13:35:12 0.0
## 10 -7139.69 1.0 2174 13:35:12 0.0
## 11 -7138.33 1.0 2174 13:35:12 0.0
## 12 -7138.00 1.0 2174 13:35:12 0.0
## 13 -7137.97 1.0 2174 13:35:12 0.0
## 14 -7137.97 1.0 2174 13:35:12 0.0
## [1] "------- Analyzing trait DON -------"
## [1] "DON has multiple years of data"
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:13 2023
## LogLik Sigma2 DF wall cpu
## 1 -6247.703 137.370 2001 13:35:13 0.0
## 2 -6221.203 126.069 2001 13:35:13 0.0
## 3 -6199.139 112.958 2001 13:35:13 0.0
## 4 -6192.242 105.175 2001 13:35:13 0.0
## 5 -6191.155 102.685 2001 13:35:13 0.0
## 6 -6191.001 102.813 2001 13:35:13 0.0
## 7 -6190.992 102.769 2001 13:35:13 0.0
## [1] "------- Analyzing trait PH -------"
## [1] "PH has a single year of data"
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:13 2023
## LogLik Sigma2 DF wall cpu
## 1 -1340.569 14.0952 649 13:35:13 0.0
## 2 -1339.579 13.7940 649 13:35:13 0.0
## 3 -1338.798 13.4275 649 13:35:13 0.0
## 4 -1338.562 13.1236 649 13:35:13 0.0
## 5 -1338.561 13.1084 649 13:35:13 0.0
BLUEs_mem<-BLUEs_mem %>%
select(GENOTYPE, FHB_ME, FDK_ME, DON_ME, HD_ME, PH_ME)
pairs.panels(BLUEs_mem[,2:ncol(BLUEs_mem)], stars = T, lm=T)
Heritabilities
#make an empty data frame for the heritability estimates
h2s<-data.frame(Trait=character(), Type=character(),
Estimation=numeric(), SE=numeric())
#calculate the harmonic means for each trait
e_traits<-list()
r_traits<-list()
for(j in traits){
a<-pheno %>%
drop_na(all_of(j))
obs_per_ind<-aggregate(a[,j]~GENOTYPE,
data=a,
length)
r<-(1/mean(1/obs_per_ind$`a[, j]`))
r_traits[j]<-r
e<-c()
for(i in unique(pheno$GENOTYPE)){
a<-pheno %>%
drop_na(j) %>%
filter(GENOTYPE==all_of(i))
a<-length(unique(a$ENV))
a<-ifelse(a==0,0,1/a)
e<-rbind(e,a)
}
e_traits[j]<-(1/mean(e))
}## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(j)
##
## # Now:
## data %>% select(all_of(j))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: There was 1 warning in `filter()`.
## ℹ In argument: `GENOTYPE == all_of(i)`.
## Caused by warning:
## ! Using `all_of()` outside of a selecting function was deprecated in tidyselect
## 1.2.0.
## ℹ See details at
## <https://tidyselect.r-lib.org/reference/faq-selection-context.html>for (i in traits){
print(paste("------- Analyzing trait", i, "-------"))
n_loc<-pheno %>%
select(ENV, YEAR, LOC, all_of(i)) %>%
drop_na() %>%
distinct(ENV)
n_loc=as.numeric(nrow(n_loc))
n_rep=n_loc*2
n_year=pheno %>%
select(ENV, YEAR, LOC, all_of(i)) %>%
drop_na() %>%
distinct(YEAR)
n_year<-as.numeric(count(n_year))
e<-as.numeric(e_traits[i])
r<-as.numeric(r_traits[i])
if(((n_year)>=2)){
print(paste(i, "has multiple years of data"))
if(i=="FDK"){
#Run Model
mlm<-asreml(fixed = pheno[,i]~1,
random = ~ GENOTYPE+
ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75,
family= asr_inverse.gaussian())
#print summary of variance comps
print(summary(mlm)$varcomp)
#pull predictions
pph2<-vpredict(mlm,
h2~V2/(V2+V3+V4))
emh2<-vpredict(mlm,
h2~V2/(V2+(V3/e)+(V4/(e*r))))
#bind predictions
a<-rbind(pph2,emh2)
a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
h2s<-rbind(h2s,a)
}else
if(i=="DON"){
#Run Model
mlm<-asreml(fixed = pheno[,i]~1,
random = ~ (GENOTYPE)+
ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75,
family= asr_inverse.gaussian())
#print summary of variance comps
print(summary(mlm)$varcomp)
#pull predictions
pph2<-vpredict(mlm,
h2~V2/(V2+V3+V4))
emh2<-vpredict(mlm,
h2~V2/(V2+(V3/e)+(V4/(e*r))))
#bind predictions
a<-rbind(pph2,emh2)
a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
h2s<-rbind(h2s,a)
}else{
#Run Model
mlm<-asreml(fixed = pheno[,i]~1,
random = ~ (GENOTYPE)+
ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75)
#print summary of variance comps
print(summary(mlm)$varcomp)
#pull predictions
pph2<-vpredict(mlm,
h2~V2/(V2+V3+V4))
emh2<-vpredict(mlm,
h2~V2/(V2+(V3/e)+(V4/(e*r))))
#bind predictions
a<-rbind(pph2,emh2)
a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
h2s<-rbind(h2s,a)
}
}else{
print(paste(i, "has a single year of data"))
mlm<-asreml(fixed = pheno[,i]~1,
random = ~ (GENOTYPE)+
ENV+
GENOTYPE:ENV,
residual = ~idv(units),
data = pheno,
maxit=75)
print(summary(mlm)$varcomp)
pph2<-vpredict(mlm,
h2~V2/(V2+V3+V4))
emh2<-vpredict(mlm,
h2~V2/(V2+(V3/e)+(V4/(e*r))))
a<-rbind(pph2,emh2)
a<-cbind(data.frame(Trait=c(i, i)), data.frame(Type=c("Per-Plot", "Entry-Mean")), a)
h2s<-rbind(h2s,a)
}
}## [1] "------- Analyzing trait HD -------"
## [1] "HD has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:25 2023
## LogLik Sigma2 DF wall cpu
## 1 -2670.166 1.0 1647 13:35:26 0.0 (1 restrained)
## Log-likelihood decreased to -2954.94 ; trying with reduced updates (0.89959)
## 2 -2607.889 1.0 1647 13:35:26 0.0
## 3 -2299.376 1.0 1647 13:35:26 0.0
## 4 -2040.728 1.0 1647 13:35:26 0.0
## 5 -2010.474 1.0 1647 13:35:26 0.0
## 6 -2002.298 1.0 1647 13:35:26 0.0
## 7 -1998.614 1.0 1647 13:35:26 0.0
## 8 -1997.200 1.0 1647 13:35:26 0.0
## 9 -1996.769 1.0 1647 13:35:26 0.0
## 10 -1996.694 1.0 1647 13:35:26 0.0
## 11 -1996.689 1.0 1647 13:35:26 0.0
## component std.error z.ratio bound %ch
## ENV 245.739732 199.1495133 1.233946 P 0.5
## GENOTYPE 3.605674 0.4143780 8.701411 P 0.0
## GENOTYPE:ENV 0.607059 0.1179896 5.145023 P 0.0
## units!units 2.626815 0.1212908 21.657163 P 0.0
## units!R 1.000000 NA NA F 0.0
## [1] "------- Analyzing trait FHB -------"
## [1] "FHB has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:26 2023
## LogLik Sigma2 DF wall cpu
## 1 -3876.512 1.0 1974 13:35:26 0.0
## 2 -3132.065 1.0 1974 13:35:26 0.0
## 3 -2339.901 1.0 1974 13:35:26 0.0
## 4 -1838.373 1.0 1974 13:35:26 0.0
## 5 -1616.057 1.0 1974 13:35:26 0.0
## 6 -1575.165 1.0 1974 13:35:26 0.0
## 7 -1572.247 1.0 1974 13:35:26 0.0
## 8 -1572.210 1.0 1974 13:35:26 0.0
## 9 -1572.210 1.0 1974 13:35:26 0.0
## component std.error z.ratio bound %ch
## ENV 0.5111042 0.36461208 1.401775 P 0
## GENOTYPE 1.9613936 0.21721428 9.029763 P 0
## GENOTYPE:ENV 0.3003627 0.05016457 5.987546 P 0
## units!units 1.1396337 0.04955938 22.995319 P 0
## units!R 1.0000000 NA NA F 0
## [1] "------- Analyzing trait FDK -------"
## [1] "FDK has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:26 2023
## LogLik Sigma2 DF wall cpu
## 1 -25196.97 1.0 2369 13:35:26 0.0
## 2 -20635.67 1.0 2369 13:35:26 0.0
## 3 -15530.46 1.0 2369 13:35:26 0.0 (1 restrained)
## 4 -11815.59 1.0 2369 13:35:26 0.0 (1 restrained)
## 5 -9558.00 1.0 2369 13:35:26 0.0 (1 restrained)
## 6 -8590.93 1.0 2369 13:35:26 0.0 (1 restrained)
## 7 -8136.15 1.0 2369 13:35:26 0.0 (1 restrained)
## 8 -7924.54 1.0 2369 13:35:26 0.0 (1 restrained)
## 9 -7835.23 1.0 2369 13:35:26 0.0
## 10 -7800.37 1.0 2369 13:35:26 0.0
## 11 -7795.71 1.0 2369 13:35:26 0.0
## 12 -7795.29 1.0 2369 13:35:26 0.0
## 13 -7795.25 1.0 2369 13:35:26 0.0
## 14 -7795.25 1.0 2369 13:35:26 0.0
## component std.error z.ratio bound %ch
## ENV 84.30977 53.697038 1.570101 P 0
## GENOTYPE 282.11502 30.766069 9.169681 P 0
## GENOTYPE:ENV 35.34972 6.574869 5.376490 P 0
## units!units 179.14713 7.110011 25.196462 P 0
## units!R 1.00000 NA NA F 0
## [1] "------- Analyzing trait DON -------"
## [1] "DON has multiple years of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:26 2023
## LogLik Sigma2 DF wall cpu
## 1 -21231.93 1.0 2196 13:35:26 0.0
## 2 -17367.32 1.0 2196 13:35:26 0.0
## 3 -13048.60 1.0 2196 13:35:26 0.0
## 4 -9918.73 1.0 2196 13:35:26 0.0
## 5 -8036.34 1.0 2196 13:35:26 0.0
## 6 -7256.48 1.0 2196 13:35:26 0.0
## 7 -6932.90 1.0 2196 13:35:26 0.0
## 8 -6810.85 1.0 2196 13:35:26 0.0
## 9 -6770.62 1.0 2196 13:35:26 0.0
## 10 -6760.85 1.0 2196 13:35:26 0.0
## 11 -6758.98 1.0 2196 13:35:26 0.0
## 12 -6758.57 1.0 2196 13:35:26 0.0
## 13 -6758.52 1.0 2196 13:35:26 0.0
## 14 -6758.52 1.0 2196 13:35:26 0.0
## component std.error z.ratio bound %ch
## ENV 77.14525 49.088534 1.571553 P 0.1
## GENOTYPE 102.31633 12.302329 8.316826 P 0.0
## GENOTYPE:ENV 52.62786 5.225355 10.071633 P 0.0
## units!units 102.48418 4.276448 23.964793 P 0.0
## units!R 1.00000 NA NA F 0.0
## [1] "------- Analyzing trait PH -------"
## [1] "PH has a single year of data"
## Model fitted using the sigma parameterization.
## ASReml 4.1.0 Sun May 7 13:35:26 2023
## LogLik Sigma2 DF wall cpu
## 1 -3067.202 1.0 842 13:35:26 0.0
## 2 -2656.660 1.0 842 13:35:26 0.0
## 3 -2215.257 1.0 842 13:35:26 0.0
## 4 -1928.060 1.0 842 13:35:26 0.0
## 5 -1793.587 1.0 842 13:35:26 0.0
## 6 -1765.862 1.0 842 13:35:26 0.0
## 7 -1763.549 1.0 842 13:35:26 0.0
## 8 -1763.517 1.0 842 13:35:26 0.0
## 9 -1763.517 1.0 842 13:35:26 0.0
## component std.error z.ratio bound %ch
## ENV 2.290131 3.3092394 0.6920414 P 0
## GENOTYPE 24.369899 3.0176516 8.0757828 P 0
## GENOTYPE:ENV 3.241033 1.0384778 3.1209455 P 0
## units!units 13.092888 0.8586791 15.2477070 P 0
## units!R 1.000000 NA NA F 0kable(h2s, row.names = F, caption="Heritability Estimates and Their Resultant Standard Errors")| Trait | Type | Estimate | SE |
|---|---|---|---|
| HD | Per-Plot | 0.5271801 | 0.0308343 |
| HD | Entry-Mean | 0.9345059 | 0.0103703 |
| FHB | Per-Plot | 0.5766447 | 0.0287133 |
| FHB | Entry-Mean | 0.9570499 | 0.0066519 |
| FDK | Per-Plot | 0.5680795 | 0.0280565 |
| FDK | Entry-Mean | 0.9690813 | 0.0049837 |
| DON | Per-Plot | 0.3974555 | 0.0304737 |
| DON | Entry-Mean | 0.9018159 | 0.0132475 |
| PH | Per-Plot | 0.5987128 | 0.0356496 |
| PH | Entry-Mean | 0.8825275 | 0.0216826 |
write.csv(h2s,
"heritabilities.csv",
row.names = FALSE)Make the rQTL file
#manipulate the genetic file to look like r qtl wants it
rownames<-colnames(geno)
geno<-as.data.frame(t(geno))
geno<-cbind(rownames,geno)
colnames(geno)<-geno[1,]
geno<-geno[-1,]
rownames(geno)<-c()
geno<-rownames_to_column(geno, var="dummy") #create a dummy column to sort later on
geno$dummy<-as.numeric(geno$dummy)
colnames(geno)[2]<-"GENOTYPE"
#merge all the blues
BLUEs_all<-left_join(BLUEs_mem, BLUEs_loc, by="GENOTYPE")
#merge the pheno and geno files
insertme<-full_join(BLUEs_all, geno, by="GENOTYPE")
insertme<-insertme[order(insertme$dummy),] #order by dummy
insertme<-insertme[1,]
pheno_geno<-merge(BLUEs_all, geno, by="GENOTYPE", all.x = T)
pheno_geno<-rbind(insertme,pheno_geno)
pheno_geno<-pheno_geno[order(pheno_geno$dummy),] #order by dummy
pheno_geno<-pheno_geno %>% filter(!is.na(dummy)==T) #get rid of any genotype that is not present
pheno_geno[1,(1:30)]<-"" #replace any NAs with blank data
dropme<-c("13955-AGS-2026","13955-JAMESTOWN", "13955-NCAG11", "13955-NC13-20076", "13955-GA06493-13LE6")
pheno_geno<-pheno_geno[!(pheno_geno$GENOTYPE %in% dropme),]
pheno_geno<-pheno_geno %>%
select(-dummy)
colnames(pheno_geno)<-gsub("FHB", "VR", colnames(pheno_geno))
#write out rqtl file
write.csv(pheno_geno,
"rqtl_input_file.csv",
row.names = F)Make a pairs plot of traits
pairs_plot_dat<-BLUEs_mem
colnames(pairs_plot_dat)<-c("Genotype",
"Visual Rating (1-9)",
"Fusarium Damaged Kernels (%)",
"Deoxynivalenol Content (PPM)",
"Heading Date (Days)",
"Plant Height (cm)")
jpeg(filename = "pairs_plot.jpg",
width = 9,
height = 9,
units = "in",
res = 320)
pairs.panels(pairs_plot_dat[,2:ncol(pairs_plot_dat)],
hist.col = "gray",
lm = TRUE,
stars = TRUE,
digits = 2,
density = FALSE,
ellipses = FALSE)
dev.off()## png
## 2Build Linkage Map
#read in file
cross_file<-read.cross("csv",
file="rqtl_input_file.csv",
genotypes=c("A","B","-"),
alleles=c("A","B"),
crosstype="dh")## --Read the following data:
## 187 individuals
## 2954 markers
## 29 phenotypes## Warning in summary.cross(cross): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
## (Perhaps it is in basepairs?)## --Cross type: dh#Pull distorted markers
gt<-geno.table(cross_file)
dropme<-rownames(gt[gt$P.value<0.0001,])
cross_file<-drop.markers(cross_file, dropme)
#remove junk
remove(dropme, gt)#create map with cM instead of BP
cross_file<- mstmap(cross_file,
pop.type="dh",
id = "GENOTYPE",
chr=c("1A",
"1B",
"1D",
"2A",
"2B",
"2D",
"3A",
"3B",
"3D",
"4A",
"4B",
"4D",
"5A",
"5B",
"5D",
"6A",
"6B",
"6D",
"7A",
"7B",
"7D"),
anchor = TRUE,
detectBadData = TRUE,
bychr = TRUE,
miss.thresh = .15,
mvest.bc = FALSE)#jitter map to improve marker distance
cross_file<-jittermap(cross_file)
summary(cross_file)## Doubled haploids
##
## No. individuals: 187
##
## No. phenotypes: 29
## Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100
## 100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9
## 100 100 100 100
##
## No. chromosomes: 85
## Autosomes: 1A.1 1A.2 1B.1 1B.2 1D.1 1D.2 1D.3 1D.4 1D.5 1D.6 2A.1
## 2A.2 2A.3 2A.4 2A.5 2B.1 2B.2 2B.3 2B.4 2D.1 2D.2 2D.3
## 2D.4 2D.5 2D.6 2D.7 2D.8 3A.1 3A.2 3A.3 3A.4 3B.1 3B.2
## 3B.3 3B.4 3B.5 3D.1 3D.2 3D.3 3D.4 3D.5 4A.1 4A.2 4B.1
## 4B.2 4B.3 4B.4 4D.1 4D.2 5A.1 5A.2 5A.3 5A.4 5B.1 5B.2
## 5D.1 5D.2 5D.3 5D.4 6A.1 6A.2 6A.3 6B.1 6B.2 6B.3 6D.1
## 6D.2 7A.1 7A.2 7A.3 7A.4 7A.5 7B.1 7B.2 7B.3 7B.4 7B.5
## 7D.1 7D.2 7D.3 7D.4 7D.5 7D.6 7D.7 7D.8
##
## Total markers: 2248
## No. markers: 111 1 133 1 57 1 1 1 1 1 133 1 1 1 1 137 4 1 1 27 1 2 1
## 1 2 1 4 177 1 1 1 218 1 1 1 1 1 14 1 1 1 120 4 70 1 1 1
## 1 12 154 1 1 1 97 5 18 1 1 4 124 2 1 123 1 1 50 2 221 1
## 1 1 1 123 1 1 1 1 40 1 3 1 1 2 3 1
## Percent genotyped: 92
## Genotypes (%): AA:54.8 BB:45.2names(cross_file$geno)## [1] "1A.1" "1A.2" "1B.1" "1B.2" "1D.1" "1D.2" "1D.3" "1D.4" "1D.5" "1D.6"
## [11] "2A.1" "2A.2" "2A.3" "2A.4" "2A.5" "2B.1" "2B.2" "2B.3" "2B.4" "2D.1"
## [21] "2D.2" "2D.3" "2D.4" "2D.5" "2D.6" "2D.7" "2D.8" "3A.1" "3A.2" "3A.3"
## [31] "3A.4" "3B.1" "3B.2" "3B.3" "3B.4" "3B.5" "3D.1" "3D.2" "3D.3" "3D.4"
## [41] "3D.5" "4A.1" "4A.2" "4B.1" "4B.2" "4B.3" "4B.4" "4D.1" "4D.2" "5A.1"
## [51] "5A.2" "5A.3" "5A.4" "5B.1" "5B.2" "5D.1" "5D.2" "5D.3" "5D.4" "6A.1"
## [61] "6A.2" "6A.3" "6B.1" "6B.2" "6B.3" "6D.1" "6D.2" "7A.1" "7A.2" "7A.3"
## [71] "7A.4" "7A.5" "7B.1" "7B.2" "7B.3" "7B.4" "7B.5" "7D.1" "7D.2" "7D.3"
## [81] "7D.4" "7D.5" "7D.6" "7D.7" "7D.8"#remove small linkage groups of a specific length
chrLengths<-c() #set up list for chromosome lengths
for (i in cross_file$geno){
chrLengths <- c(chrLengths, length(i$map))
}
#create data frame with number of markers per linkage group
chrLengths<-as.data.frame(cbind(names(cross_file$geno), chrLengths))
colnames(chrLengths)<-c("chr", "n_markers")
chrLengths$n_markers<-as.numeric(chrLengths$n_markers)
chrLengths<-filter(chrLengths, n_markers > 5)
#new number of linkage groups
nrow(chrLengths)## [1] 21chrLengths## chr n_markers
## 1 1A.1 111
## 2 1B.1 133
## 3 1D.1 57
## 4 2A.1 133
## 5 2B.1 137
## 6 2D.1 27
## 7 3A.1 177
## 8 3B.1 218
## 9 3D.2 14
## 10 4A.1 120
## 11 4B.1 70
## 12 4D.2 12
## 13 5A.1 154
## 14 5B.1 97
## 15 5D.1 18
## 16 6A.1 124
## 17 6B.1 123
## 18 6D.1 50
## 19 7A.1 221
## 20 7B.1 123
## 21 7D.1 40#subset the map with the linkage groups that pass the filter
cross_file<-subset(cross_file, chr = chrLengths$chr)
cross_file<-jittermap(cross_file)
#show summary and write file with output
summary(cross_file)## Doubled haploids
##
## No. individuals: 187
##
## No. phenotypes: 29
## Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100
## 100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9
## 100 100 100 100
##
## No. chromosomes: 21
## Autosomes: 1A.1 1B.1 1D.1 2A.1 2B.1 2D.1 3A.1 3B.1 3D.2 4A.1 4B.1
## 4D.2 5A.1 5B.1 5D.1 6A.1 6B.1 6D.1 7A.1 7B.1 7D.1
##
## Total markers: 2159
## No. markers: 111 133 57 133 137 27 177 218 14 120 70 12 154 97 18 124
## 123 50 221 123 40
## Percent genotyped: 92
## Genotypes (%): AA:54.7 BB:45.3write.cross(cross_file, "cross_file_without_parents",format = "csv")
#remove junk
remove(i, chrLengths)
#check map metrics
names(cross_file$pheno)## [1] "GENOTYPE" "VR_ME" "FDK_ME" "DON_ME" "HD_ME"
## [6] "PH_ME" "VR_KIN2019" "FDK_KIN2019" "DON_KIN2019" "HD_KIN2020"
## [11] "VR_KIN2020" "FDK_KIN2020" "DON_KIN2020" "PH_KIN2020" "HD_LW2019"
## [16] "VR_LW2019" "FDK_LW2019" "DON_LW2019" "HD_LW2020" "VR_LW2020"
## [21] "FDK_LW2020" "DON_LW2020" "PH_LW2020" "FDK_VIR2019" "DON_VIR2019"
## [26] "HD_VIR2020" "VR_VIR2020" "FDK_VIR2020" "DON_VIR2020"names(cross_file)## [1] "geno" "pheno"names(cross_file$geno)## [1] "1A.1" "1B.1" "1D.1" "2A.1" "2B.1" "2D.1" "3A.1" "3B.1" "3D.2" "4A.1"
## [11] "4B.1" "4D.2" "5A.1" "5B.1" "5D.1" "6A.1" "6B.1" "6D.1" "7A.1" "7B.1"
## [21] "7D.1"nphe(cross_file)## [1] 29nind(cross_file)## [1] 187nchr(cross_file)## [1] 21summary(cross_file)## Doubled haploids
##
## No. individuals: 187
##
## No. phenotypes: 29
## Percent phenotyped: 100 100 100 100 100 100 97.9 97.9 97.9 100 100 100 100
## 100 98.9 98.9 98.9 98.9 100 100 100 100 100 98.9 98.9
## 100 100 100 100
##
## No. chromosomes: 21
## Autosomes: 1A.1 1B.1 1D.1 2A.1 2B.1 2D.1 3A.1 3B.1 3D.2 4A.1 4B.1
## 4D.2 5A.1 5B.1 5D.1 6A.1 6B.1 6D.1 7A.1 7B.1 7D.1
##
## Total markers: 2159
## No. markers: 111 133 57 133 137 27 177 218 14 120 70 12 154 97 18 124
## 123 50 221 123 40
## Percent genotyped: 92
## Genotypes (%): AA:54.7 BB:45.3#rename linkage groups
names(cross_file$geno)<-c("1A",
"1B",
"1D",
"2A",
"2B",
"2D",
"3A",
"3B",
"3D",
"4A",
"4B",
"4D",
"5A",
"5B",
"5D",
"6A",
"6B",
"6D",
"7A",
"7B",
"7D")
#plot images of map
plotMap(cross_file, horizontal = FALSE, shift = FALSE, main="NC13-20076 Double Haploid Linkage Map")
geno.image(cross_file, main = "")
Check Map
#check orentation of markers
mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)## $Markers
## [1] "character"
##
## $Chr
## [1] "character"
##
## $BP
## [1] "character"
##
## $cM
## [1] "numeric"pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
geom_point()+
facet_wrap(~Chr,
ncol = 3,
scales = "free")+
labs(title = "Centimorgan Position vs. Megabase Pair",
x = "Centimorgan (cM) Position",
y = "Megabase Pair (Mbp) Position")
for (i in names(cross_file$geno)){
q<-subset(pos, pos$Chr==i)
if (q[1,3]>q[(nrow(q)),3]){
cross_file<-flip.order(cross_file, i)
print(paste("Chromsome",i,"has been flipped"))
remove(q)
}else{
print(paste("Chromsome",i,"is in correct order"))
remove(q)
}
}## [1] "Chromsome 1A is in correct order"
## [1] "Chromsome 1B is in correct order"
## [1] "Chromsome 1D is in correct order"
## [1] "Chromsome 2A is in correct order"
## [1] "Chromsome 2B is in correct order"
## [1] "Chromsome 2D is in correct order"
## [1] "Chromsome 3A is in correct order"
## [1] "Chromsome 3B is in correct order"
## [1] "Chromsome 3D is in correct order"
## [1] "Chromsome 4A is in correct order"
## [1] "Chromsome 4B is in correct order"
## [1] "Chromsome 4D is in correct order"
## [1] "Chromsome 5A is in correct order"
## [1] "Chromsome 5B is in correct order"
## [1] "Chromsome 5D is in correct order"
## [1] "Chromsome 6A is in correct order"
## [1] "Chromsome 6B is in correct order"
## [1] "Chromsome 6D is in correct order"
## [1] "Chromsome 7A is in correct order"
## [1] "Chromsome 7B is in correct order"
## [1] "Chromsome 7D has been flipped"mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)## $Markers
## [1] "character"
##
## $Chr
## [1] "character"
##
## $BP
## [1] "character"
##
## $cM
## [1] "numeric"pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
geom_point()+
facet_wrap(~Chr,
ncol = 3,
scales = "free")+
labs(title = "Centimorgan Position vs. Megabase Pair",
x = "Centimorgan (cM) Position",
y = "Megabase Pair (Mbp) Position")
#remove bad markers
removeme<-c("S1A_527395239", "S1A_208718176", "S7A_621573544")
cross_file<-drop.markers(cross_file, removeme)
#check one last time
mn<-markernames(cross_file)
pos<-find.markerpos(cross_file,mn)
pos<-rownames_to_column(pos,var="BP_Location")
pos<-pos %>% separate(BP_Location, into =c("a","b"), sep = "_")
colnames(pos)<-c("trash","BP","Chr","cM")
pos<-pos %>% select(Chr, BP, cM)
pos<-cbind(mn,pos)
colnames(pos)<-c("Markers", "Chr", "BP", "cM")
lapply(pos, class)## $Markers
## [1] "character"
##
## $Chr
## [1] "character"
##
## $BP
## [1] "character"
##
## $cM
## [1] "numeric"pos[,3:4]<-lapply(pos[,3:4], as.numeric)
ggplot(data = pos, aes(x=cM, y=BP/1000000))+
geom_point()+
facet_wrap(~Chr,
ncol = 3,
scales = "free")+
labs(title = "Centimorgan Position vs. Megabase Pair",
x = "Centimorgan (cM) Position",
y = "Megabase Pair (Mbp) Position")
ggsave(plot = last_plot(),
filename = "cM_vs_bp_comparison.jpg",
units = "in",
width = 8,
height = 12,
dpi = 320)
#remove
remove(mn, pos, removeme)#calculate marker probabilities
cross_file<-calc.genoprob(cross_file,
step=2.0,
off.end=0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#calculate marker probabilities for a simulated genotype
cross_file<-sim.geno(cross_file,
n.draws = 128,
step = 2,
off.end = 0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#detect the number of cores of the computer processor
ncores<-as.numeric(detectCores())
#set the number of permutations
nperms=1000
#### Run Heading Date ####
traits<-names(cross_file$pheno)[c(grep("HD", names(cross_file$pheno)),grep("PH", names(cross_file$pheno)))]
#make results vector
results<-list()
#run initial interval mapping and pull out QTL
for (i in traits){
#announce
print(paste("------------ Interval Mapping of", i,"------------"))
#perform IM with multiple imputation method
print("Interval mapping...")
scans<-scanone(cross_file,
pheno.col = i,
addcovar = NULL,
model = "normal",
method = "hk")
print("Done")
#perform IM permutations mapping to define significance threshold
print("Permutational interval mapping...")
perms<-scanone(cross_file,
pheno.col = i,
addcovar = NULL,
model = "normal",
method = "hk",
n.perm = nperms,
n.cluster = ncores-1) #set threshold
print("Done")
#plot QTL Scan
print("Plotting...")
threshold<-summary(perms, alpha=0.05)
plot(scans,main=paste("IM for", i))
abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
#print plot
jpeg(paste("Scan_IM_", i, ".jpg", sep = ""),
width = 11,
height = 4,
units = "in",
res = 320)
threshold<-summary(perms, alpha=0.05)
plot(scans,main=paste("IM for", i))
abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
dev.off()
print("Done")
#show the peak markers for QTL
print("Defining QTL...")
qtl<-summary(scans, perm=perms,
lodcolum=1, alpha=0.05)
#rename QTL to identify which scan they came from
qtl$name=paste("IM_", qtl$chr,"-pos-",qtl$pos,sep="")
print("Done")
#Define the QTL locations and effects
print("Drawing QTL...")
colnames(scans)<-c("chr","pos","lod")
#place outputs in lists
results$IM$scan[[i]]<-scans #place the scan a the list
results$IM$perms[[i]]<-perms #place the permutations a the list
results$IM$threshold[[i]]<-threshold #place the thresholds a the list
#set up objects for defining QTL
c<-qtl[,1] #define the chromosomes where QTL are found
p<-qtl[,2] #define the positions of the QTL
a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
#Make new cross with genome subset
a<-sim.geno(a,
n.draws = 128,
step = 2,
off.end = 0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#make a QTL object from that subset
madeqtl<-makeqtl(a,
c,
p,
qtl.name = qtl[,4],
what = c("prob"))
#place that QTL object in a list
results$IM$qtl[[i]]<-madeqtl
print("Done")
#announce
print("Running drop-one QTL analysis to check significance...")
#test the significance of those QTL using drop one analysis
qtlfit<-fitqtl(cross_file,
qtl=results$IM$qtl[[i]],
pheno.col = i,
model="normal",
method = "hk")
#remove insignificant qtl
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
madeqtl<-dropfromqtl(madeqtl,
qtl.name = a$name)
results$IM$qtl[[i]]<-madeqtl
remove(a,c,p,scans,perms,threshold)
}else{
remove(a,c,p,scans,perms,threshold)
}
#set vectors outside the loop
j=1
#make initial check object
qtl_check<-results$IM$qtl[[i]]
#run MQM until no significant peaks!
repeat{
#announce
print(paste("-------- Performing Additional QTL Scan for Trait", i,"--------"))
mqm<-addqtl(cross = cross_file,
pheno.col = i,
qtl=qtl_check,
covar = NULL,
method = "hk")
results[[paste("MQM", j, sep="")]]$scan[[i]]<-mqm
#plot QTL Scan
print("Plotting...")
plot(mqm,main=paste("MQM", j, "for", i))
abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
print("Done")
#write out picture to directory
jpeg(paste("Scan_",paste("MQM", j, sep=""),"_", i, ".jpg", sep = ""),
width = 11,
height = 4,
units = "in",
res = 320)
plot(mqm,main=paste("MQM", j, "for", i))
abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
dev.off()
#make a qtl object
qtl<-summary(mqm,
perm=results$IM$perms[[i]],
lodcolum=1,
alpha=0.05)
if(nrow(qtl)==0){
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
print(paste("There were no new QTL identified for", i, "aborting loop"))
break
}else{
#rename QTL to identify which scan they came from
qtl$name=paste(paste("MQM", j, "_", sep=""), qtl$chr,"-pos-",qtl$pos,sep="")
#Defining QTL
print("Drawing new QTL")
c<-qtl[,1] #define the chromosomes where QTL are found
p<-qtl[,2] #define the positions of the QTL
a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
#simulate the genome for that subset
a<-sim.geno(a,
n.draws = 128,
step = 2,
off.end = 0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#make a QTL object from that subset
madeqtl<-makeqtl(a,
c,
p,
qtl.name = qtl[,4],
what = c("prob"))
#place that QTL object in a list
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
print("Done")
#announce
print("Running drop-one QTL analysis...")
#test the significance of those QTL using dropone analysis
qtlfit<-fitqtl(cross_file,
qtl=madeqtl,
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
if(!is.null(qtlfit$result.drop)){
print("Checking drop-one analysis for significance")
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
a<-tidyr::separate(data = a,
col = "name",
into = c("chr", "trash", "pos"),
sep = "-")
a$pos<-as.numeric(a$pos)
a$chr<-gsub(paste("MQM", j, "_", sep=""),"",a$chr)
madeqtl<-dropfromqtl(madeqtl,
chr = a$chr,
pos = a$pos)
madeqtl<-addtoqtl(cross_file,
qtl = madeqtl,
chr = qtl_check$chr,
pos = qtl_check$pos,
qtl.name = qtl_check$name)
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
remove(a,c,p)
}else{
print(paste("There was only one new QTL identified for ",
i,
" in ",
paste("MQM", j, sep=""),
"... Checking significance!",
sep = ""))
if(qtlfit$result.full[1,6]<0.05){
print("New QTL is significant, placing in total model!")
madeqtl<-addtoqtl(cross_file,
qtl = madeqtl,
chr = qtl_check$chr,
pos = qtl_check$pos,
qtl.name = qtl_check$name)
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
}else{
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
print(paste("There were no new QTL identified for", i, "aborting loop!"))
break
}
}
}
qtl_check<-results[[paste("MQM", j, sep="")]]$qtl[[i]]
j=j+1
}
#check
if(j-1==0){
#pull and make final qtl object
finalqtl<-results$IM$qtl[[i]]
}else{
#pull and make final qtl object
finalqtl<-results[[paste("MQM", j-1, sep="")]]$qtl[[i]]
}
#check for significance
qtlfit<-fitqtl(cross_file,
qtl=finalqtl,
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
#rename the qtl
a<-finalqtl
a<-data.frame(summary(a))
a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
finalqtl$name=a$name
remove(a)
#print
print("Filtering insignificant markers...")
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
finalqtl<-dropfromqtl(finalqtl,
qtl.name = a$name)
results$final$qtl_unrefine[[i]]<-finalqtl
remove(a, finalqtl)
}else{
remove(a, finalqtl)
}
#print
print("Refining QTL position...")
#refine
results$final$qtl_refine[[i]]<-refineqtl(cross = cross_file,
pheno.col = i,
qtl = results$final$qtl_unrefine[[i]],
verbose = FALSE,
method = "hk")
#rename the qtl
a<-results$final$qtl_refine[[i]]
a<-data.frame(summary(a))
a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
results$final$qtl_refine[[i]]$name=a$name
remove(a)
#check for significance
qtlfit<-fitqtl(cross_file,
qtl=results$final$qtl_refine[[i]],
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
finalqtl<-results$final$qtl_refine[[i]]
#print
print("Filtering insignificant markers...")
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
finalqtl<-dropfromqtl(finalqtl,
qtl.name = a$name)
results$final$qtl_refine_filt[[i]]<-finalqtl
remove(a, finalqtl)
}else{
remove(a, finalqtl)
}
droplodint<-c()
#create QTL 1.5 drop LOD intervals
for(z in results$final$qtl_refine_filt[[i]]$name){
a<-data.frame(name=z)
a<-tidyr::separate(data = a,
col = "name",
into = c("scan", "other"),
sep = "_",
remove = FALSE)
a<-tidyr::separate(data = a,
col = "other",
into = c("chr", "trash", "pos"),
sep = "-",
remove = TRUE)
b<-results[[a$scan]]$scan[[i]]
c<-find.pseudomarker(cross = cross_file,
chr = a$chr,
pos = as.numeric(a$pos))
d<-lodint(results = b,
chr = a$chr,
qtl.index = c,
drop = 1.5,
lodcolumn = 1,
expandtomarkers = TRUE)
d<-data.frame(name=z,
trait=i,
chr=unique(d$chr),
start_marker=rownames(d)[1],
start=d[1,2],
mid_marker=rownames(d)[2],
mid_pos=d[2,2],
stop_marker=rownames(d)[3],
stop=d[3,2])
droplodint<-rbind(droplodint, d)
remove(a,b,c,d)
}
#place in results
results$final$qtl_support_intervals[[i]]<-droplodint
#remove
remove(droplodint)
}## [1] "------------ Interval Mapping of HD_ME ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_ME --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_ME --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for HD_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_KIN2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_KIN2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_KIN2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_KIN2020 --------"
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for HD_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_LW2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = NULL, : Dropping 2 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "There was only one new QTL identified for HD_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = NULL, : Dropping 2 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for HD_LW2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_LW2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_LW2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for HD_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of HD_VIR2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait HD_VIR2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for HD_VIR2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait HD_VIR2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for HD_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_ME ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_ME --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for PH_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait PH_ME --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for PH_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_KIN2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait PH_KIN2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for PH_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of PH_LW2020 ------------"
## [1] "Interval mapping..."
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."
## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait PH_LW2020 --------"
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for PH_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."#rename results
results_HD_PH<-resultsPull Marker Covariates
covars<-rbind(results_HD_PH$final$qtl_support_intervals$HD_ME,
results_HD_PH$final$qtl_support_intervals$PH_ME)
covar_markers<-data.frame(GENOTYPE=cross_file$pheno$GENOTYPE)
for(i in 1:nrow(covars)){
a<-covars[i,]
b<-fill.geno(cross = cross_file,
method = c("no_dbl_XO"),
map.function = "kosambi")
b<-pull.geno(cross = b,
chr = a$chr)
a<-find.marker(cross = cross_file,
chr = a$chr,
pos = a$mid_pos)
b<-as.data.frame(b[,a])
colnames(b)=a
covar_markers<-cbind(covar_markers, b)
remove(a,b)
}Visualize HD and PH QTL
#make object for ggplot
traits<-names(cross_file$pheno)[c(grep("HD", names(cross_file$pheno)),grep("PH", names(cross_file$pheno)))]
#make object
hd_ph_ggplot<-c()
for(i in traits){
a<-results_HD_PH$final$qtl_support_intervals[[i]]
hd_ph_ggplot<-rbind(hd_ph_ggplot, a)
remove(a)
}
d<-results_HD_PH$IM$scan$HD_ME %>%
filter(chr %in% hd_ph_ggplot$chr)
hd_ph_ggplot<-hd_ph_ggplot %>%
mutate(chr=paste(chr, trait, sep = "_"))
library(ggrepel)
library(patchwork)
finalplot<-list()
j=1
for(i in unique(d$chr)){
a<-hd_ph_ggplot[grep(i, hd_ph_ggplot$chr),]
a$trait<-factor(a$trait, levels = unique(hd_ph_ggplot$trait))
b<-d[d$chr==i,]
b<-b[-grep("loc", rownames(b)),]
b<-rownames_to_column(b, var = "marker")
Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
if(j==1){
c<-ggplot(data=b, aes(x=chr, y=pos))+
geom_point(size=2)+
geom_line(linewidth=1)+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("HD_ME" = "blue",
"HD_KIN2020" = "dodgerblue",
"HD_LW2019" = "darkslategray",
"HD_LW2020" = "cyan",
"HD_VIR2020" = 'cadetblue',
"PH_ME" = 'violet',
"PH_KIN2020" = 'darkorchid',
"PH_LW2020"="mediumorchid"),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())+
labs(title = i,
color = "Trait Within Environment",
y = "cM")
if(length(Label)>0){
c<-c+
geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
}
}else{
c<-ggplot(data = b, aes(x=chr, y=pos))+
geom_point(size=2)+
geom_line(linewidth=1)+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("HD_ME" = "blue",
"HD_KIN2020" = "dodgerblue",
"HD_LW2019" = "darkslategray",
"HD_LW2020" = "cyan",
"HD_VIR2020" = 'cadetblue',
"PH_ME" = 'violet',
"PH_KIN2020" = 'darkorchid',
"PH_LW2020"="mediumorchid"),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank())+
labs(title = i,
color = "Trait Within Environment",
y = "cM")
if(length(Label)>0){
c<-c+
geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
}
}
finalplot[[paste("chr",i,sep="_")]]<-c
remove(a,b,c)
j=j+1
}
a<-(finalplot$chr_1B+finalplot$chr_2A+finalplot$chr_3A+
finalplot$chr_4A+finalplot$chr_5A+finalplot$chr_6A+
finalplot$chr_7A+finalplot$chr_7B+finalplot$chr_7D)+
plot_layout(guides = "collect", nrow = 1)+
plot_annotation(title = "1.5-LOD Support Intervals",
subtitle = "Plant Height (PH) and Heading Date (HD)")
a
ggsave(plot = a,
filename = "support_intervals_hd_and_ph.jpg",
width = 14,
height = 8,
dpi = 320,
units = "in")Visualize HD and PH QTL for ME only
#make object for ggplot
traits<-c("HD_ME", "PH_ME")
#make object
hd_ph_ggplot<-c()
for(i in traits){
a<-results_HD_PH$final$qtl_support_intervals[[i]]
hd_ph_ggplot<-rbind(hd_ph_ggplot, a)
remove(a)
}
d<-results_HD_PH$IM$scan$HD_ME %>%
filter(chr %in% hd_ph_ggplot$chr)
hd_ph_ggplot<-hd_ph_ggplot %>%
mutate(trait=ifelse(trait=="HD_ME", "Heading Date",
ifelse(trait=="PH_ME", "Plant Height", "ERROR"))) %>%
mutate(chr=paste(chr, trait, sep = "_"))
library(ggrepel)
library(patchwork)
finalplot<-list()
j=1
for(i in unique(d$chr)){
a<-hd_ph_ggplot[grep(i, hd_ph_ggplot$chr),]
a$trait<-factor(a$trait, levels = unique(hd_ph_ggplot$trait))
b<-d[d$chr==i,]
b<-b[-grep("loc", rownames(b)),]
b<-rownames_to_column(b, var = "marker")
Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
if(j==1){
c<-ggplot(data=b, aes(x=chr, y=pos))+
geom_point(size=2)+
geom_line(linewidth=1)+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("Heading Date" = "blue",
"Plant Height" = 'violet'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())+
labs(title = i,
color = "Trait",
y = "cM")
if(length(Label)>0){
c<-c+
geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
}
}else{
c<-ggplot(data = b, aes(x=chr, y=pos))+
geom_point(size=2)+
geom_line(linewidth=1)+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("Heading Date" = "blue",
"Plant Height" = 'violet'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank())+
labs(title = i,
color = "Trait",
y = "cM")
if(length(Label)>0){
c<-c+
geom_point(data = d[rownames(d) %in% Label,], aes(x=chr, y=pos), color="green", size=4, shape=17)
}
}
finalplot[[paste("chr",i,sep="_")]]<-c
remove(a,b,c)
j=j+1
}
a<-(finalplot$chr_2A+finalplot$chr_4A+finalplot$chr_5A+
finalplot$chr_6A+finalplot$chr_7B+finalplot$chr_7D)+
plot_layout(guides = "collect", nrow = 1)+
plot_annotation(title = "1.5-LOD Support Intervals",
subtitle = "Plant Height (PH) and Heading Date (HD)")
a
ggsave(plot = a,
filename = "support_intervals_hd_and_ph_ME_only.jpg",
width = 14,
height = 8,
dpi = 320,
units = "in")Run Mapping for Disease Traits
#set the number of permutations
nperms=1000
#### Run disease traits ####
traits<-names(cross_file$pheno)[c(grep("VR", names(cross_file$pheno)),grep("FDK", names(cross_file$pheno)),grep("DON", names(cross_file$pheno)))]
#make results vector
results<-list()
#run initial interval mapping and pull out QTL
for (i in traits){
#announce
print(paste("------------ Interval Mapping of", i,"------------"))
#perform IM with multiple imputation method
print("Interval mapping...")
scans<-scanone(cross_file,
pheno.col = i,
addcovar = covar_markers[,2:ncol(covar_markers)],
model = "normal",
method = "hk")
print("Done")
#perform IM permutations mapping to define significance threshold
print("Permutational interval mapping...")
perms<-scanone(cross_file,
pheno.col = i,
addcovar = covar_markers[,2:ncol(covar_markers)],
model = "normal",
method = "hk",
n.perm = nperms,
n.cluster = ncores-1) #set threshold
print("Done")
#plot QTL Scan
print("Plotting...")
threshold<-summary(perms, alpha=0.05)
plot(scans,main=paste("IM for", i))
abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
#print plot
jpeg(paste("Scan_IM_", i, ".jpg", sep = ""),
width = 11,
height = 4,
units = "in",
res = 320)
threshold<-summary(perms, alpha=0.05)
plot(scans,main=paste("IM for", i))
abline(h=threshold, lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
dev.off()
print("Done")
#show the peak markers for QTL
print("Defining QTL...")
qtl<-summary(scans, perm=perms,
lodcolum=1, alpha=0.05)
#rename QTL to identify which scan they came from
qtl$name=paste("IM_", qtl$chr,"-pos-",qtl$pos,sep="")
print("Done")
#Define the QTL locations and effects
print("Drawing QTL...")
colnames(scans)<-c("chr","pos","lod")
#place outputs in lists
results$IM$scan[[i]]<-scans #place the scan a the list
results$IM$perms[[i]]<-perms #place the permutations a the list
results$IM$threshold[[i]]<-threshold #place the thresholds a the list
#set up objects for defining QTL
c<-qtl[,1] #define the chromosomes where QTL are found
p<-qtl[,2] #define the positions of the QTL
a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
#Make new cross with genome subset
a<-sim.geno(a,
n.draws = 128,
step = 2,
off.end = 0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#make a QTL object from that subset
madeqtl<-makeqtl(a,
c,
p,
qtl.name = qtl[,4],
what = c("prob"))
#place that QTL object in a list
results$IM$qtl[[i]]<-madeqtl
print("Done")
#announce
print("Running drop-one QTL analysis to check significance...")
#test the significance of those QTL using drop one analysis
qtlfit<-fitqtl(cross_file,
qtl=results$IM$qtl[[i]],
pheno.col = i,
model="normal",
method = "hk")
#remove insignificant qtl
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
madeqtl<-dropfromqtl(madeqtl,
qtl.name = a$name)
results$IM$qtl[[i]]<-madeqtl
remove(a,c,p,scans,perms,threshold)
}else{
remove(a,c,p,scans,perms,threshold)
}
#set vectors outside the loop
j=1
#make initial check object
qtl_check<-results$IM$qtl[[i]]
#run MQM until no significant peaks!
repeat{
#announce
print(paste("-------- Performing Additional QTL Scan for Trait", i,"--------"))
mqm<-addqtl(cross = cross_file,
pheno.col = i,
qtl=qtl_check,
covar = covar_markers[,2:ncol(covar_markers)],
method = "hk")
results[[paste("MQM", j, sep="")]]$scan[[i]]<-mqm
#plot QTL Scan
print("Plotting...")
plot(mqm,main=paste("MQM", j, "for", i))
abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
print("Done")
#write out picture to directory
jpeg(paste("Scan_",paste("MQM", j, sep=""),"_", i, ".jpg", sep = ""),
width = 11,
height = 4,
units = "in",
res = 320)
plot(mqm,main=paste("MQM", j, "for", i))
abline(h=results$IM$threshold[[i]], lty="dotted", lwd=1, col="#cc0000")
legend("topleft",legend=c("p=0.05"),col = c("#cc0000"),lty = "dotted")
dev.off()
#make a qtl object
qtl<-summary(mqm,
perm=results$IM$perms[[i]],
lodcolum=1,
alpha=0.05)
if(nrow(qtl)==0){
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
print(paste("There were no new QTL identified for", i, "aborting loop"))
break
}else{
#rename QTL to identify which scan they came from
qtl$name=paste(paste("MQM", j, "_", sep=""), qtl$chr,"-pos-",qtl$pos,sep="")
#Defining QTL
print("Drawing new QTL")
c<-qtl[,1] #define the chromosomes where QTL are found
p<-qtl[,2] #define the positions of the QTL
a<-subset(cross_file, chr=c) #subset the chromosomes where QTL are found
#simulate the genome for that subset
a<-sim.geno(a,
n.draws = 128,
step = 2,
off.end = 0.0,
error.prob=1.0e-4,
map.function="kosambi",
stepwidth="fixed")
#make a QTL object from that subset
madeqtl<-makeqtl(a,
c,
p,
qtl.name = qtl[,4],
what = c("prob"))
#place that QTL object in a list
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
print("Done")
#announce
print("Running drop-one QTL analysis...")
#test the significance of those QTL using dropone analysis
qtlfit<-fitqtl(cross_file,
qtl=madeqtl,
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
if(!is.null(qtlfit$result.drop)){
print("Checking drop-one analysis for significance")
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
a<-tidyr::separate(data = a,
col = "name",
into = c("chr", "trash", "pos"),
sep = "-")
a$pos<-as.numeric(a$pos)
a$chr<-gsub(paste("MQM", j, "_", sep=""),"",a$chr)
madeqtl<-dropfromqtl(madeqtl,
chr = a$chr,
pos = a$pos)
madeqtl<-addtoqtl(cross_file,
qtl = madeqtl,
chr = qtl_check$chr,
pos = qtl_check$pos,
qtl.name = qtl_check$name)
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
remove(a,c,p)
}else{
print(paste("There was only one new QTL identified for ",
i,
" in ",
paste("MQM", j, sep=""),
"... Checking significance!",
sep = ""))
if(qtlfit$result.full[1,6]<0.05){
print("New QTL is significant, placing in total model!")
madeqtl<-addtoqtl(cross_file,
qtl = madeqtl,
chr = qtl_check$chr,
pos = qtl_check$pos,
qtl.name = qtl_check$name)
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-madeqtl
}else{
results[[paste("MQM", j, sep="")]]$qtl[[i]]<-NULL
print(paste("There were no new QTL identified for", i, "aborting loop!"))
break
}
}
}
qtl_check<-results[[paste("MQM", j, sep="")]]$qtl[[i]]
j=j+1
}
#check
if(j-1==0){
#pull and make final qtl object
finalqtl<-results$IM$qtl[[i]]
}else{
#pull and make final qtl object
finalqtl<-results[[paste("MQM", j-1, sep="")]]$qtl[[i]]
}
#check for significance
qtlfit<-fitqtl(cross_file,
qtl=finalqtl,
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
#rename the qtl
a<-finalqtl
a<-data.frame(summary(a))
a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
finalqtl$name=a$name
remove(a)
#print
print("Filtering insignificant markers...")
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
finalqtl<-dropfromqtl(finalqtl,
qtl.name = a$name)
results$final$qtl_unrefine[[i]]<-finalqtl
remove(a, finalqtl)
}else{
remove(a, finalqtl)
}
#print
print("Refining QTL position...")
#refine
results$final$qtl_refine[[i]]<-refineqtl(cross = cross_file,
pheno.col = i,
qtl = results$final$qtl_unrefine[[i]],
verbose = FALSE,
method = "hk")
#rename the qtl
a<-results$final$qtl_refine[[i]]
a<-data.frame(summary(a))
a<-tidyr::separate(a, col = "name", into=c("model", "trash"), sep = "_")
a$name<-paste(a$model, "_", a$chr, "-pos-", a$pos, sep = "")
results$final$qtl_refine[[i]]$name=a$name
remove(a)
#check for significance
qtlfit<-fitqtl(cross_file,
qtl=results$final$qtl_refine[[i]],
pheno.col = i,
get.ests = T,
model="normal",
method = "hk")
finalqtl<-results$final$qtl_refine[[i]]
#print
print("Filtering insignificant markers...")
if(!is.null(qtlfit$result.drop)){
a<-as.data.frame(qtlfit$result.drop)
a$sig<-ifelse(a$`Pvalue(F)`<0.05, 1, 0)
a<-data.frame(name=rownames(a[a$sig==0,]))
finalqtl<-dropfromqtl(finalqtl,
qtl.name = a$name)
results$final$qtl_refine_filt[[i]]<-finalqtl
remove(a, finalqtl)
}else{
remove(a, finalqtl)
}
droplodint<-c()
#create QTL 1.5 drop LOD intervals
for(z in results$final$qtl_refine_filt[[i]]$name){
a<-data.frame(name=z)
a<-tidyr::separate(data = a,
col = "name",
into = c("scan", "other"),
sep = "_",
remove = FALSE)
a<-tidyr::separate(data = a,
col = "other",
into = c("chr", "trash", "pos"),
sep = "-",
remove = TRUE)
b<-results[[a$scan]]$scan[[i]]
c<-find.pseudomarker(cross = cross_file,
chr = a$chr,
pos = as.numeric(a$pos))
d<-lodint(results = b,
chr = a$chr,
qtl.index = c,
drop = 1.5,
lodcolumn = 1,
expandtomarkers = TRUE)
d<-data.frame(name=z,
trait=i,
chr=unique(d$chr),
start_marker=rownames(d)[1],
start=d[1,2],
mid_marker=rownames(d)[2],
mid_pos=d[2,2],
stop_marker=rownames(d)[3],
stop=d[3,2])
droplodint<-rbind(droplodint, d)
remove(a,b,c,d)
}
#place in results
results$final$qtl_support_intervals[[i]]<-droplodint
#remove
remove(droplodint)
}## [1] "------------ Interval Mapping of VR_ME ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait VR_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_KIN2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait VR_KIN2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_KIN2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_KIN2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_KIN2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_LW2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "There was only one new QTL identified for VR_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_LW2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of VR_LW2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for VR_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for VR_LW2020 in MQM2... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait VR_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."## Warning: Expected 3 pieces. Additional pieces discarded in 1 rows [1].## [1] "------------ Interval Mapping of VR_VIR2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait VR_VIR2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for VR_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_ME ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_KIN2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Checking drop-one analysis for significance"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_KIN2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_KIN2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_KIN2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_LW2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "There was only one new QTL identified for FDK_LW2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_LW2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_LW2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for FDK_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_VIR2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_VIR2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of FDK_VIR2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for FDK_VIR2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait FDK_VIR2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for FDK_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_ME ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for DON_ME in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_ME --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_ME aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_KIN2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 4 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait DON_KIN2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 14 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_KIN2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 4 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_KIN2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_KIN2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_KIN2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_LW2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_LW2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_LW2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."
## [1] "There was only one new QTL identified for DON_LW2020 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_LW2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_LW2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_VIR2019 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 2 individuals with missing phenotypes.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "Drawing new QTL"
## [1] "Done"
## [1] "Running drop-one QTL analysis..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "There was only one new QTL identified for DON_VIR2019 in MQM1... Checking significance!"
## [1] "New QTL is significant, placing in total model!"
## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2019 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 12 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"## [1] "There were no new QTL identified for DON_VIR2019 aborting loop"## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 2 individuals with missing phenotypes.## [1] "Filtering insignificant markers..."
## [1] "------------ Interval Mapping of DON_VIR2020 ------------"
## [1] "Interval mapping..."## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 10 individuals with missing covariates.
## [1] "Done"
## [1] "Permutational interval mapping..."
## -Running permutations via a cluster of 11 nodes.
## [1] "Done"
## [1] "Plotting..."## [1] "Done"
## [1] "Defining QTL..."
## [1] "Done"
## [1] "Drawing QTL..."
## [1] "Done"
## [1] "Running drop-one QTL analysis to check significance..."
## [1] "-------- Performing Additional QTL Scan for Trait DON_VIR2020 --------"## Warning in addqtl(cross = cross_file, pheno.col = i, qtl = qtl_check, covar = covar_markers[, : Dropping 10 individuals with missing phenotypes.
## [1] "Plotting..."## [1] "Done"
## [1] "There were no new QTL identified for DON_VIR2020 aborting loop"
## [1] "Filtering insignificant markers..."
## [1] "Refining QTL position..."
## [1] "Filtering insignificant markers..."#rename results
results_FHB_FDK_DON<-resultsVisualize FHB, FDK, and DON QTL
#make object for ggplot
traits<-names(cross_file$pheno)[c(grep("VR", names(cross_file$pheno)),grep("FDK", names(cross_file$pheno)),grep("DON", names(cross_file$pheno)))]
#make object
fhb_fdk_don_ggplot<-c()
for(i in traits){
a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
remove(a)
}
d<-results_FHB_FDK_DON$IM$scan$VR_ME %>%
filter(chr %in% as.character(unique(fhb_fdk_don_ggplot$chr)))
fhb_fdk_don_ggplot<-fhb_fdk_don_ggplot %>%
mutate(trait=gsub("FHB", "VR", trait)) %>%
mutate(chr=paste(chr, trait, sep = "_"))
library(ggrepel)
library(patchwork)
finalplot<-list()
j=1
for(i in unique(d$chr)){
a<-fhb_fdk_don_ggplot[grep(i, fhb_fdk_don_ggplot$chr),]
a$trait<-factor(a$trait, levels = unique(fhb_fdk_don_ggplot$trait))
b<-d[d$chr==i,]
b<-b[-grep("loc", rownames(b)),]
b<-rownames_to_column(b, var = "marker")
Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
if(j==1 | j==11){
c<-ggplot(data=b, aes(x=chr, y=pos))+
geom_point(size=1)+
geom_line()+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("VR_ME" = "red",
"VR_KIN2019" = "firebrick",
"VR_KIN2020" = "indianred",
"VR_LW2019" = "brown",
"VR_LW2020" = 'orangered',
"VR_VIR2020" = 'tomato',
"FDK_ME" = 'yellow',
"FDK_KIN2019" = 'darkgoldenrod',
"FDK_KIN2020" = 'gold',
"FDK_LW2019" = 'goldenrod',
"FDK_LW2020" = 'lightgoldenrod',
"FDK_VIR2019" = 'darkgoldenrod2',
"FDK_VIR2020" = 'khaki4',
"DON_ME" = 'green',
"DON_KIN2019" = 'chartreuse3',
"DON_KIN2020" = 'darkolivegreen',
"DON_LW2019" = 'darkseagreen',
"DON_LW2020" = 'olivedrab',
"DON_VIR2019" = 'palegreen',
"DON_VIR2020" = 'springgreen1'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())+
labs(title = i,
color = "Trait Within Environment",
y = "cM")
}else{
c<-ggplot(data = b, aes(x=chr, y=pos))+
geom_point(size=1)+
geom_line()+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("VR_ME" = "red",
"VR_KIN2019" = "firebrick",
"VR_KIN2020" = "indianred",
"VR_LW2019" = "brown",
"VR_LW2020" = 'orangered',
"VR_VIR2020" = 'tomato',
"FDK_ME" = 'yellow',
"FDK_KIN2019" = 'darkgoldenrod',
"FDK_KIN2020" = 'gold',
"FDK_LW2019" = 'goldenrod',
"FDK_LW2020" = 'lightgoldenrod',
"FDK_VIR2019" = 'darkgoldenrod2',
"FDK_VIR2020" = 'khaki4',
"DON_ME" = 'green',
"DON_KIN2019" = 'chartreuse3',
"DON_KIN2020" = 'darkolivegreen',
"DON_LW2019" = 'darkseagreen',
"DON_LW2020" = 'olivedrab',
"DON_VIR2019" = 'palegreen',
"DON_VIR2020" = 'springgreen1'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank())+
labs(title = i,
color = "Trait Within Environment",
y = "cM")
}
finalplot[[paste("chr",i,sep="_")]]<-c
remove(a,b,c)
j=j+1
}
a<-(finalplot$chr_1A+finalplot$chr_1B+finalplot$chr_1D+
finalplot$chr_2A+finalplot$chr_2B+finalplot$chr_2D+
finalplot$chr_3A+finalplot$chr_3B+finalplot$chr_3D+
finalplot$chr_4A+finalplot$chr_4B+finalplot$chr_4D+
finalplot$chr_5A+finalplot$chr_5B+finalplot$chr_5D+
finalplot$chr_6A+finalplot$chr_6B+
finalplot$chr_7A+finalplot$chr_7D)+
plot_layout(guides = "collect",
nrow = 2,
ncol = 10)+
plot_annotation(title = "1.5-LOD Support Intervals",
subtitle = "Visual Ratings (VR), Fusarium Damaged Kernels (FDK), and Deoxynivalenol Content (DON)")
a
ggsave(plot = a,
filename = "support_intervals_vr_fdk_and_don.jpg",
width = 20,
height = 8,
dpi = 320,
units = "in")Redraw with just ME
#make object for ggplot
traits<-c("VR_ME", "FDK_ME", "DON_ME")
#make object
fhb_fdk_don_ggplot<-c()
for(i in traits){
a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
remove(a)
}
d<-results_FHB_FDK_DON$IM$scan$VR_ME %>%
filter(chr %in% as.character(unique(fhb_fdk_don_ggplot$chr)))
fhb_fdk_don_ggplot<-fhb_fdk_don_ggplot %>%
mutate(trait=ifelse(trait=="VR_ME", "Visual Rating",
ifelse(trait=="FDK_ME", "Fusarium Damaged Kernels",
ifelse(trait=="DON_ME", "Deoxynivalenol Content", "ERROR")))) %>%
mutate(chr=paste(chr, trait, sep = "_"))
library(ggrepel)
library(patchwork)
finalplot<-list()
j=1
for(i in c("2A", "3B", "4A", "4B", "6A", "7A")){
a<-fhb_fdk_don_ggplot[grep(i, fhb_fdk_don_ggplot$chr),]
a$trait<-factor(a$trait, levels = unique(fhb_fdk_don_ggplot$trait))
b<-d[d$chr==i,]
b<-b[-grep("loc", rownames(b)),]
b<-rownames_to_column(b, var = "marker")
Label<-colnames(covar_markers)[grep(i, colnames(covar_markers))]
if(j==1){
c<-ggplot(data=b, aes(x=chr, y=pos))+
geom_point(size=1)+
geom_line()+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("Visual Rating" = "red",
"Fusarium Damaged Kernels" = 'yellow',
"Deoxynivalenol Content" = 'green'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank())+
labs(title = i,
color = "Trait",
y = "cM")
}else{
c<-ggplot(data = b, aes(x=chr, y=pos))+
geom_point(size=1)+
geom_line()+
scale_y_reverse(limits = c(max(d$pos),0))+
scale_color_manual(values = c("Visual Rating" = "red",
"Fusarium Damaged Kernels" = 'yellow',
"Deoxynivalenol Content" = 'green'),
drop=FALSE)+
geom_segment(data = a, aes(x=chr, xend=chr, y=start, yend=stop, color=trait), linewidth=2)+
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.y = element_blank())+
labs(title = i,
color = "Trait",
y = "cM")
}
finalplot[[paste("chr",i,sep="_")]]<-c
remove(a,b,c)
j=j+1
}
a<-(finalplot$chr_2A+finalplot$chr_3B+finalplot$chr_4A+finalplot$chr_4B+finalplot$chr_6A+finalplot$chr_7A)+
plot_layout(guides = "collect",
nrow = 1)+
plot_annotation(title = "1.5-LOD Support Intervals",
subtitle = "Visual Ratings (VR), Fusarium Damaged Kernels (FDK), and Deoxynivalenol Content (DON)")
a
ggsave(plot = a,
filename = "support_intervals_vr_fdk_and_don_ME_only.jpg",
width = 10,
height = 8,
dpi = 320,
units = "in")Percent variance in ME
fhb_fdk_don_ggplot<-c()
for(i in traits){
a<-results_FHB_FDK_DON$final$qtl_support_intervals[[i]]
fhb_fdk_don_ggplot<-rbind(fhb_fdk_don_ggplot, a)
remove(a)
}
qtl_info<-c()
for(i in c("2A", "3B", "4A", "4B", "6A", "7A")){
a<-fhb_fdk_don_ggplot[fhb_fdk_don_ggplot$chr==i,]
vr_effects<-fitqtl(cross = cross_file,
pheno.col = "VR_ME",
covar = covar_markers[,2:ncol(covar_markers)],
qtl = results_FHB_FDK_DON$final$qtl_refine_filt$VR_ME,
get.ests = TRUE)
fdk_effects<-fitqtl(cross = cross_file,
pheno.col = "FDK_ME",
covar = covar_markers[,2:ncol(covar_markers)],
qtl = results_FHB_FDK_DON$final$qtl_refine_filt$FDK_ME,
get.ests = TRUE)
don_effects<-fitqtl(cross = cross_file,
pheno.col = "DON_ME",
covar = covar_markers[,2:ncol(covar_markers)],
qtl = results_FHB_FDK_DON$final$qtl_refine_filt$DON_ME,
get.ests = TRUE)
b<-as.data.frame(vr_effects$result.drop)[rownames(vr_effects$result.drop) %in% a$name,]
b<-rownames_to_column(b, var = "name")
b<-b %>% left_join(a, by="name") %>% filter(trait=="VR_ME")
c<-as.data.frame(fdk_effects$result.drop)[rownames(fdk_effects$result.drop) %in% a$name,]
c<-rownames_to_column(c, var = "name")
c<-c %>% left_join(a, by="name") %>% filter(trait=="FDK_ME")
d<-as.data.frame(don_effects$result.drop)[rownames(don_effects$result.drop) %in% a$name,]
d<-rownames_to_column(d, var = "name")
d<-d %>% left_join(a, by="name") %>% filter(trait=="DON_ME")
e<-as.data.frame(summary(vr_effects)$est)[rownames(summary(vr_effects)$est) %in% a$name,]
e<-rownames_to_column(e, var = "name")
b<-b %>% left_join(e, by="name") %>% filter(trait=="VR_ME")
f<-as.data.frame(summary(fdk_effects)$est)[rownames(summary(fdk_effects)$est) %in% a$name,]
f<-rownames_to_column(f, var = "name")
c<-c %>% left_join(f, by="name") %>% filter(trait=="FDK_ME")
g<-as.data.frame(summary(don_effects)$est)[rownames(summary(don_effects)$est) %in% a$name,]
g<-rownames_to_column(g, var = "name")
d<-d %>% left_join(g, by="name") %>% filter(trait=="DON_ME")
h<-rbind(b,c,d)
qtl_info<-rbind(qtl_info, h)
remove(a,b,c,d,e,f,g,h,vr_effects,fdk_effects,don_effects)
}## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "VR_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "FDK_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.## Warning in fitqtl(cross = cross_file, pheno.col = "DON_ME", covar =
## covar_markers[, : The qtl object doesn't contain imputations; using
## method="hk".## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 10 individuals with missing phenotypes.qtl_info<-qtl_info %>%
arrange(chr, mid_pos, trait) %>%
separate(start_marker, into=c("junk1", "start_bp"), sep = "_", remove = FALSE) %>%
separate(stop_marker, into=c("junk2", "stop_bp"), sep = "_", remove = FALSE) %>%
mutate(`Proximal Mbp Position`= round(as.numeric(start_bp)/1000000, 2),
`Distal Mbp Position` = round(as.numeric(stop_bp)/1000000, 2),
newname = paste("QFhb.nc-", chr, sep = ""),
`%var`=`%var`/100) %>%
select(newname,
trait,
chr,
start_marker,
start,
`Proximal Mbp Position`,
stop_marker,
stop,
`Distal Mbp Position`,
est,
SE,
`%var`) %>%
rename(QTL=newname,
Trait=trait,
Chromosome=chr,
`Proximal Flanking Marker`=start_marker,
`Proximal cM Position`=start,
`Distal Flanking Marker`=stop_marker,
`Distal cM Position`= stop,
`Estimated Effect`=est,
`Standard Error`=SE,
`Percent Variance`=`%var`)
write.csv(qtl_info,
"final_qtl_information.csv",
row.names = FALSE)